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ABSTRACT 

A numerical analysis of combustion instabilities that induce flashback in a lean, premixed, 
prevaporized dump combustor is performed. KTVA-II, a finite volume CFD code for the 
modeling of transient, multidimensional, chemically reactive flows, serves as the principal 
analytical tool. The experiment of Proctor and Tien [13] is used as a reference for 
developing the computational model. An experimentally derived combustion instability 
mechanism is presented on the basis of the observations of Proctor and Tien and other 
investigators of instabilities in low-speed (M<0.1) dump combustors. The analysis comprises 
two independent procedures that begin from a calculated stable flame: the first is a linear 
increase of the equivalence ratio and the second is the linear decrease of the inflow velocity. 
The objective is to observe changes in the aerothermochemical features of the flow field 
prior to flashback. It is found that only the linear increase of the equivalence ratio elicits 
a calculated flashback result. Though this result does not exhibit large-scale coherent 
vortices in the turbulent shear layer coincident with a flame flickering mode as observed 
experimentally, there are interesting acoustic effects which are resolved quite well in the 
calculation. A discussion of the k - e turbulence model used by KtVA-II is prompted by the 
absence of combustion instabilities in the model as the inflow velocity is linearly decreased. 
Finally, recommendations are made for further numerical analysis that may improve 
correlation with experimentally observed combustion instabilities. 
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* OVERVIEW 

A dump combustor is characterized by 1) the sudden expansion of a premixed, 
prevaporized fuel/air mixture over a rearward facing step or dump plane, and 2) combustion 
of this mixture across a flame front that is seated on the edge of the step. Combustion 
instabilities which disturb the flame have been investigated in many experiments using low- 
speed (M<0.1) axisymmetric or rectangular dump combustors. One such experiment, 
conducted at the NASA Lewis Research Center by Proctor and T’ien [13], is used as a basis 
in this computational analysis of combustion instabilities leading to flashback, a condition 
in which the flame migrates upstream into the premixing region. In addition to expounding 
on this statement of purpose. Chapter 1 examines a few experiments which both predate and 
postdate that by Proctor and T’ien. A combustion instability mechanism evolving from 
these and other experimental investigations is outlined at the end of the chapter. It is 
important to understand each step in this mechanism since it is referred to repeatedly. 

Chapter 2 is devoted to describing the analytical methods that are employed in this study. 
It begins by presenting results of the Proctor and T’ien experiment that are used in 
developing the numerical model. This is followed by a description of the ICEd-ALE finite 
volume method implemented in the KIVA-II computer program. A summary of the 
governing equations of motion, turbulence, and chemical reaction solved by KIVA-II is also 
presented. Before setting up a computational model of the Proctor and T’ien experiment, 
two preliminary tests of the KIVA-II code are performed. A brief assessment of the results 
of these tests is made at the end of the chapter with reference to Figures 2 and 3. 

Chapter 3 describes how the computational model is set up, beginning with the generation 
of the mesh based on the geometry of the Proctor and T’ien experimental rig. Details of the 
boundary conditions invoked at the inflow and outflow planes are presented. Also described 
is the law-of-the-wall boundary layer prescription that is used to model the velocity profile 
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and heat flux at the solid boundaries. Chapter 3 closes with an identification of the initial 
constituents and thermodynamic state of the fluid in the mesh. 

Chapter 4 begins with a description of the cold flow calculation. Then the difficulties that 
arise during the calculation of ignition are discussed with comparisons made to the ignition 
procedure implemented by Proctor and T’ien. Following this is a description of the method 
used to calibrate the chemical kinetics model in KIVA-II so that flashback is calculated at 
approximately the same equivalence ratio as observed by Proctor and T’ien. Finally, the 
two-part procedure used to investigate combustion instabilities prior to flashback is 
described. 

The results of this two-part procedure are presented and discussed in Chapter 5. The first 
part, which entails the linear increase of the equivalence ratio 1 from a calculated stable 
flame, does yield a flashback result. However, the second part of the procedure, which is the 
linear decrease of the inflow velocity from the same stable flame calculation, does not 
induce combustion instability or flashback in the numerical model. Suggested reasons for 
the disparity in these results are discussed. 

Chapter 6 presents a summary of the experimentally observed reactive flow phenomena that 
the KIVA-1T model simulates well. Also presented are the disparities between the calculated 
combustion instabilities and the experimentally observed mechanism outlined at the end of 
Chapter 1 . Finally, recommendations are made for improving the acoustic and chemical 
modeling capabilities of the KF/A-II program. These recommendations would increase 
data storage and computing time, but should result in improved agreement with 
experimental observations. 


1. In this analysis, the equivalence ratio is defined on a mass basis as the ratio of the actual to 
stoichiometric fuel/dry-air ratios. 
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. CHAPTER 1 
INTRODUCTION 


1.1 Statement of Purpose 

Oxides of nitrogen (NO x ) are known to be detrimental to the ozone layer. Stringent 
emission standards for these compounds have resulted in the investigation of new gas 
turbine combustor designs by NASA and private aeronautics companies. One such design 
is a lean, premixed, prevaporized (LPP) dump combustor. The development of this 
particular gas turbine combustor has been hindered by the occurrence of combustion 
instabilities which can cause the flame front to migrate into the premixing region, a 
phenomenon called flashback. Three redeeming features of the LPP dump combustor are 
simplicity of design, high combustion efficiency, and reduced NO x emissions. 

In addition to its aeronautical application, O.I. Smith, et.al. [1] illustrate by experiment the 
potential of the LPP dump combustor for waste incineration. Long residence times of high 
temperature combustion products in the recirculation zone formed behind the dump plane 
offer a means for ‘high destruction and removal efficiencies’ of hazardous waste. 
Understanding the instability mechanisms leading to flashback is of particular importance 
in this application since stable, coherent recirculation zones are required to ensure complete 
incineration. 

Culick [2] notes that the term ‘combustion instability’ commonly used in the literature is 
misleading. He reasons that combustion processes themselves, though sometimes 
oscillatory, are in most instances stable. But when they are coupled with convective 
oscillations that are acoustic, hydrodynamic, or thermodynamic in origin, sustained 
pressure fluctuations can occur in the flow field, resulting in noise, vibration, and flame 
oscillation which are characterized as combustion instabilities. 
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Three mechanisms are repeatedly discussed in the literature as drivers of combustion 
instability in an LPP dump combustor: 

(1) acoustic disturbances that result from the fact that variable speed flow of an 
inviscid gas through a dump combustor may excite acoustic modes, each 
with their own frequency; 

(2) hydrodynamic instabilities ascribed to the growth, coupling, and convection 
of large-scale, coherent vortices in the shear layer downstream of the step; 

(3) unsteady heat release oscillations resulting from combustion wave 
propagation within the cores of large-scale vortices. 

The objective of this study is to investigate the exact nature of the flow field pressure 
oscUlations and other aerothermochemical conditions existing prior to and during flashback 
for low speed flows (M<0.1) in a simple LPP dump combustor using KIVA-II, a recently 
developed computational fluid dynamics code [3], An improved understanding of these 
flow field features may lead to clearer insight into which of the three instability 
mechanisms, or combination thereof, are principal drivers under given operating 
conditions. 

1.2 Literature Review 

Many experimental studies of combustion instabilities in an LPP dump combustor have 
been published. Some focus on the high speed (M>1) instabilities that impede ramjet 
performance. However, for this analysis only those experiments that investigate instability 
in the combustion of low speed flows (M<0.1 and 10 3 <Re<10 4 ) in an LPP dump combustor 
are considered. 

Improved understanding of combustion instability has been derived largely from recent 
insight on the vortical dynamics of nonreacting turbulent shear layers. In the mid- 1970 s, 
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several investigators identified large-scale 1 , coherent vortices as principal mechanisms in 
the development and growth of turbulent shear layers [4-6]. Later, Ganji and Sawyer [7] 
were among the first to study the development of large-scale vortices in both non-reacting 
and reacting turbulent flows in a dump combustor. By varying the diameter and location of 
a trip wire in the laminar boundary layer of the premixing region, they observe that the 
formation and growth rate of coherent vortices in the shear layer downstream is impaired 
when transition to turbulence occurs in the premixer boundary layer. In addition, they note 
that there is little difference in the incipient development and growth rate of vortices in 
reacting and non-reacting shear layers. Finally, they observe that flame structure and 
propagation is influenced by developing vortices since they entrain reactants at the top of 
the shear layer and induce combustion through subsequent mixing with hot products. 

Keller, et.al. [8] investigate large-scale vortex movements and growth rates as instability 
mechanisms leading to flashback. They distinguish three levels of instability as the fuel 
mass flow rate is increased to the lean flashback limit. Each is associated with a pressure 
profile of distinct amplitude and frequency which is evidenced by noise characterized as 
humming, buzzing, and chucking. Accompanying each sound change is an increased level 
of flame oscillation laterally across the combustion zone. The humming is ascribed to 
longitudinal standing waves associated with the natural frequencies of the dump combustor, 
an acoustic effect. The buzzing and chucking are more energetic and are said to be produced 
by the growth and convection of large-scale vortices in the shear layer, a hydrodynamic 
effect. When chucking occurs, the convective oscillations produced by the action of large- 
scale vortices in the initial region of the shear layer induce flow reversals in the boundary 
layer just upstream of the step. This enables the flame front to detach from the comer of the 
step and migrate upstream through the premixer boundary layer, giving rise to flashback. 

1 ‘Large-scale’ connotes vortical diameters that are many times larger than the length scale of turbulent 
eddies. These vortices span an identifiable, and often major, portion of the width of the turbulent shear 
layer. 
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The experimental results of Smith and Zukoski [9] reveal that large-scale vortices form only 
at resonant frequencies of the combustor and that combustion within these vortices as they 
grow and convect downstream induces a periodic heat release rate. They note that, in 
accordance with the Rayleigh criterion [10], when a constructive phase relationship exists 
between flow field pressure oscillations and the heat release rate, the acoustic field is 
excited. If enough energy is added to the acoustic field, large velocity fluctuations result 
which drive the flame off the comer of the step. This is most likely analogous to the 
‘chucking’ effect in the Keller, et.al. experiment. 

Schadow, et.al. [11] describe a driving mechanism for combustion instability that is similar 
to that of Smith and Zukoski. In addition, they propose that combustion occurs first at the 
circumference of any given large-scale vortex and propagates inward, resulting in the 
periodic heat release. Schadow, et.al. also note that a turbulent shear layer has associated 
with it a range of hydrodynamic instability frequencies. In its initial region, the shear layer 
develops instability waves whose frequencies are a function of the flow structure in the 
boundary layer upstream of the step as described in references [4-6] and [7]. Further along 
in the shear layer are instabilities that are subharmonics of the initial shear layer instability 
frequency. The size that vortices achieve as they convect downstream is determined by the 
frequency at which the most constructive phase relationship occurs between the shear layer 
instabilities and the system resonance mode. Schadow, et.al. conclude that decoupling these 
oscillatory mechanisms, particularly when the development of large-scale vortices is 
favored, could be a means of damping the large pressure oscillations that drive combustion 
instability. (It is important to note that, in contrast to the other experiments described herein, 
Schadow, et.al. investigate a diffusion flame in a dump combustor. However, in considering 
the qualitative effects of large-scale vortices on combustion instability, it is assumed that 
there is little or no distinction between premixed and diffusion flames.) 
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Logan, et.al. [12] have developed a very interesting experimental dump combustor in which 
they are able to adjust the length (l c ) to height (h c ) aspect ratio of the combustion chamber 
using movable ceramic plugs at the exit (Figure la). They distinguish three modes of 
instability at various combinations of combustion chamber geometry and equivalence ratio 
(Figure lb). The first is a low-frequency (30-50Hz) ‘chugging’ mode in which a single 
vortex pair periodically develops at the dump plane and grows to fill the entire combustor 
cavity, precluding the establishment of a stable recirculation zone. The other two modes, 
distinguished as ‘A’ and ‘B’, are at much higher frequencies (mode A:600-700Hz, mode 
B:350-600Hz) and elicit the simultaneous occurrence of a few or several vortex pairs in the 
combustion chamber. A coherent recirculation zone is maintained at these high frequency 
modes. Between the shaded regions in Figure lb are quiescent conditions in which the 
recirculation zone is unperturbed. One of the most revealing insights of this experimental 
study is the effect of the acoustic field on the appearance of combustion instabilities and on 
the potential for modulating them with adjustments to the combustion chamber geometry. 

An experimental investigation of flame flashback in an LPP dump combustor is also done 
by Proctor and T’ien [13] at NASA’s Lewis Research Center. They propose a flashback 
mechanism similar to that of Keller, et.al. which entails the propagation of the flame 
upstream through the premixer boundary layer. At first, this appears to be a reaffirmation 
of the ‘classical’ flashback mechanism of Lewis and von Elbe which relates the flame speed 
to the mean velocity in the boundary layer of a Bunsen burner tube [21]. However, Proctor 
and T’ien observe that lower equivalence ratios are required for incipient flashback as the 
premixer velocity is increased. Since flame speed increases with equivalence ratio, Proctor 
and T’ien conclude that a different flashback mechanism exists in a dump combustor. 
Consistent with their observations, they propose that flow reversals in the premixer 
boundary layer elicited by high amplitude combustor pressure oscillations actuate upstream 
flame propagation. This experiment is discussed in more detail in the next chapter. 
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1.3 Combustion Instability Mechanism Revealed by Experiment 

The increased amplitude of combustor pressure oscillations is consistently observed to be 
coincident with flame flickering and, in the most extreme case, flashback in a dump 
combustor. Some of the experiments described in the previous section have examined the 
vortical dynamics of reacting turbulent shear layers in conjunction with system resonance 
modes to explain the development of these high amplitude pressure oscillations. The 
following outline serves as a summarized combustion instability mechanism based on the 
experimental work cited in the previous section and in references [18-20]: 

(1) an acoustic field exists with fundamental and harmonic frequencies unique 
to the combustor geometry; 

(2) turbulent shear layer instability waves exist whose frequencies are a function 
of the velocity profile and flow structure of the boundary layer just upstream 
of the dump plane; 

(3) a constructive phase relationship arises at one or more frequencies common 
to the shear layer and acoustic field, leading to the development of coherent 
vortices; 

(4) combustion within the vortices as they grow and convect downstream results 
in a periodic heat release rate which intensifies the acoustic field; 

(5) the stronger acoustic field favors development of larger vortices which, in 
turn, amplify the heat release oscillations; 

(6) self-excited oscillations in the ‘closed loop’ mechanism of steps (l)-(5) 
cause high amplitude pressure oscillations that ultimately drive the flame 
upstream into the premixer with its leading edge in the premixer boundary 
layer. 

An understanding of the parameters controlling the occurrence of step (3) is essential since 
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large-scale vortices elicit oscillations in the heat release rate. Experiment has cited fuel 
type, equivalence ratio, premixer Reynolds number, and combustor geometry as controlling 
parameters in the development of coherent vortices in the shear layer. In this analysis, 
combustor geometry and fuel type are fixed while the equivalence ratio and the premixer 
Reynolds number are varied independently from a stable flame solution. It is important that 
each step in the experimentally observed instability mechanism outlined above be 
understood since it is frequently referenced in this analysis. 
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. CHAPTER 2 
ANALYTICAL METHODS 


2.1 Experimental Basis 

The experimental work by Proctor and T’ien serves as a basis for this computational 
analysis. In this experiment, inlet air temperature, premixer wall temperature, and average 
premixer velocity are varied parametrically to determine their individual effects on the 
minimum equivalence ratio required for sustained flashback. The results of their parametric 
study are: 

(1) the effect of inlet air temperature is ‘slight’; 

(2) premixer wall temperature is said not to be an important parameter in the 
occurrence of flashback; 

(3) as the average premixer velocity is increased, lower equivalence ratios are 
required for flashback. 

Proctor and T’ien conclude that flashback is induced by flow reversals in the boundary layer 
at the premixer exit. These flow reversals are most likely generated by pressure oscillations 
in the combustion chamber which are observed to reach a maximum amplitude just before 
flashback. Flame flickering prior to flashback, perhaps caused by large-scale vortices 
shedding from the step, is considered a visible sign of the increased amplitude of the 
oscillations (Figure 5a+b). 

2.2 Computational Methods 

The KIVA-II computer program developed at Los Alamos National Laboratory is used. It 
solves the governing equations for transient, multidimensional, chemically reactive flow 
using a combined Implicit Continuous-flow Eulerian and Arbitrary Lagrangian-Eulerian 
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(ICEd-ALE) finite volume method. The synergistic pairing of the ICE and ALE methods 
enable KIVA-II to model a wide range of flows at all speeds (ICE) in arbitrary geometries 
(ALE). By itself, the ALE method has two advantages: 1) a deformable mesh comprised of 
hexahedrons whose vertices are arbitrarily specified functions of time, and 2) quick 
convergence of the explicit solution procedures for high speed flows. A disadvantage, 
however, is that for low speed flows the Courant sound speed stability criterion is extremely 
restrictive on the time step. By combining the implicit calculation of pressure in the ICE 
method with the ALE method, the stringency of the Courant condition is removed so that 
low speed flows can be calculated efficiently. 

The ICEd-ALE method comprises three phases. Phase A is an explicit calculation that 
solves the Lagrangian conservation equations for mass and momentum as a couple. The 
Lagrangian volumes that result are used to solve the energy equation, yielding the time n+l 
internal energies. Finally in Phase A, the time n+l pressures are computed from the 
equation of state. For high speed flows in which a completely explicit calculation is 
acceptable, Phase A would, by itself, give a solution. For low speed flows, however, the 
implicit calculation of Phase B is necessary so that pressure signals may traverse more than 
one cell in a time step. Using the Lagrangian velocities determined in Phase A as a first 
guess. Phase B invokes mass conservation as it implicitly updates the pressure, density, and 
Lagrangian velocities according to the ICE method. Then, just as in Phase A, time n+l 
internal energy and pressure are obtained from the solution of the energy and state 
equations. 

Since Phases A and B actuate a purely Lagrangian description of the flow field, mesh 
vertices move at the local fluid velocity, preventing advection across cell boundaries. In 
Phase C, called the rezone phase, explicit solution methods are used to fix the flow field, 
move the vertices to new locations, and then remap or ‘rezone’ the flow field onto the new 
mesh by advecting mass, momentum, and energy across cell boundaries. In this way Phase 
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C effectively undoes some or all of the changes to the shape of the mesh that occur in the 
Lagrangian calculations of Phases A and B. Since a continuous flow device with 
immovable boundaries is considered in this analysis, Phase C maintains a fixed 
computational mesh by moving vertices back to their original positions every time step for 
a purely Eulerian description of the flow field. 

In the ICEd-ALE method, equations are differenced in integral form using the volume of a 
typical computational cell as a control volume. Divergence terms are transformed to surface 
integrals using the divergence theorem [3]. 


2.3 Governing Equations 


An overview of the governing equations solved by the ICEd-ALE finite volume scheme 
follows. Note that V is the three dimensional del operator defined as: 


)d r ; d 


v = *£ + >£ +w - 


3 _ 

dz 


Conservation Equations 
Continuity: 

!^ + ^.(p m 5) = v.} m + p£, 

for the mth species, where p m is the mass density of species m, u is the fluid velocity, J m 
is the diffusive mass flux of species m, and p c m is the rate of change of p m due to chemical 
reaction. The diffusive mass flux is given by Fick’s Law: 

K = povy m 

where p is the total mass density and Y m is the species mass fraction. The mass diffusivity 
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D is assumed to be the same for all species since turbulence is the dominant mechanism for 
mass diffusion. It is computed using the Schmidt number Sc and viscosity p: 


D = 


pSc * 


By summing over all species, the total mass density equation is: 

^ + V. (pfi) = 0 

since mass is conserved in diffusion and chemical reaction. 


Momentum: 

+ V • (pw«) = -Vp-Vpk + V • a 

where p is the fluid pressure, and k is the specific turbulent kinetic energy. The viscous 
stress tensor is Newtonian and is expressed as: 

g = p[Vu+ (Vm) 7 ] +W*uU 

where p and \ are the shear and bulk viscosity coefficients, U is the unit dyadic tensor, and 

superscript T denotes the transpose. The viscosity coefficients are given by: 

^ k 1 
M- i i air +C |iP £ 

and 

i 2 

where e is the specific turbulent kinetic energy dissipation rate and is an empirical 
turbulence constant equalling 0.09. A Sutherland formula is used for the molecular 
viscosity of air with A \ and A 2 constant: 

_ A,r 3/2 
“ t + a 2 * 
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Internal Energy Transport 


The transport of internal energy p/ is governed by: 


3(p/) 

dt 


+ V»(p«7) = -pV • u - V • q + pe+ Q c 


where £) c is a source term due to chemical heat release. The heat flux q, is gained by 
summing the effects of heat conduction and enthalpy transport by mass diffusion: 

m 

where K is the thermal conductivity, T is the fluid temperature, and h m is the species 
specific enthalpy. The thermal conductivity is determined using the Prandtl number Pr : 

k ~T7 

where c p is the specific heat at constant pressure. Finally, pe is the rate of internal energy 
production from turbulent kinetic energy dissipation. 

Turbulence Equations 

A two-equation k-z turbulence model is used. It solves the following pair of equations for 
the turbulent kinetic energy pfc, and the turbulent kinetic energy dissipation rate pe, 
respectively: 


^ + V • (puk) = kV • u + g:Vu + V • -pe 

at 3 L rr * J 


9pe 

IT 


+ V«(p uz) = [-c|^ £i -c e3 } ] p eV . + V • + ^[c £ a:Vu -c^pe] 


where c , c , c , Pr t , and Pr art empirical constants having the following values [14] 

V e 2 E 3 * t 
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c f = 1.44 

^ i 


Pr k = 1.0 


c f = 1.92 
fc 2 


Pr e = 1.3 


Cp = - 1.0 . 


With values for fc and e, the turbulent length scale L is computed directly from: 


,, 3/2 


L = c 


e 


where 




f >r t( c e 2 c e,) 


1/2 


Chemical Equations 

Gaseous mixtures are assumed to behave as ideal gases and therefore conform to the 
equation of state: 

p = 

V m m / 

where W is the molecular weight of species m and R u is the universal gas constant. 


$ 


i 


’•5 


■■*> E- 

r 1 


Chemical kinetics model: 


Using propane (C 3 Hg) as fuel, four kinetic reactions are considered: 


C 3 tf 8 + 5 0 2 3 C0 2 + 4tf 2 0 


0 2 + 2N 2 *+2N+2N0 


20 2 + N 2 <r>20 + 2N0 


N 2 + 20H-*2H + 2N0 


A global reaction rate mechanism is employed to approximate the effects of the many 
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elementary steps associated with each of these reactions. The global rate expression is: 

■ • 

where a' mr and b ^ may be the stoichiometric coefficients of reactants m and products n in 
the four kinetic reactions considered or may be empirical constants determined from 
experiment or detailed kinetics calculations [15]. The forward and backward reaction rate 
coefficients *y r and k br for reaction r are determined from the Arrhenius expression: 



The temperature exponent the activation energy E a , and the pre-exponential term A, are 
also determined from correlation with experimental results in reference [15]. 

Equilibrium chemistry model: 

Six equilibrium reactions are considered: 

H 2 <-> 2H 

0 2 <r> 20 

N 2 <r*2N 

0 2 + H 2 <-> 20H 
0 2 + 2H 2 0<^40H 
0 2 + 2CO 2C0 2 • 

The rates of each equilibrium reaction above are determined implicitly from: 



where v ' and v" are the stoichiometric coefficients of species m on the left and right 
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side, respectively, of each reaction and K[ (T) is the known equilibrium constant. 

2.4 Preliminary Test of Computational Model 

KIVA-II was originally developed to model the confined flows of internal combustion 
engines. Subsequent code modifications have added an inflow/outflow capability for gas 
turbine applications. To substantiate these changes, preliminary analyses of two flow 
regimes are conducted. The first is a non-reacting turbulent flow over a step. Using identical 
Reynolds numbers, turbulence parameters, and computational meshes, the results of the 
KIVA-II program are qualitatively similar to those of the FIDAP code developed and 
marketed by Fluid Dynamics International, Inc. (Figure 2a+b). It should be noted that the 
disparity in the numerical values accompanying the KIVA-II and FIDAP graphical output 
for this calculation occurs because KIVA-II uses the CGS system of units while FIDAP uses 
nondimensionalized variables. A direct quantitative comparison is difficult since the non- 
dimensionalization parameters used in the FIDAP simulation are not specified. However, it 
can be seen that the KIVA-II and FIDAP simulations yield the same order of magnitude for 
the difference between maximum and minimum values of corresponding flow field 
quantities. 

The Keller, et.al. experiment, introduced in section 1.2, is used as a second test. Figure 3a 
shows cinematographic schlieren records of a stable flame at an equivalence ratio of 0.57. 
With corresponding mesh geometry, boundary and initial conditions, fluid properties, 
Reynolds number, and turbulence parameters, KIVA-II yields the solution shown in Figure 
3b. This solution is at steady state with virtually no combustor pressure oscillations, which 
contrasts with the highly transient nature of the experimental flame. The k - e turbulence 
model used in KIVA-II does not appear to simulate the chaotic undulation of the flame in 
the latter part of the combustor. This shortcoming is discussed further in Chapter 5 with 
regard to the combustion model of the Proctor and T’ien experiment. 
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CHAPTER 3 
MODEL DEFINITION 


3.1 Computational Mesh 

Calculation is made in two dimensions on the centrally located x-z plane (y=0.0in.) of the 
Proctor and T’ien experimental dump combustor (Figure 4a). Flow in this plane is 
symmetric about x=0.0in. and, for this analysis, calculation occurs only in the lower half 
(Figure 4b). To properly model the incompressible flow resistance, the cross-sectional area 
ratio of the combustor to exhaust aperture in the computational mesh (A*/B*) closely 
approximates the corresponding ratio for the experimental dump combustor (A/B): 

£ = 5.093 £ = 5.333 . 

B B 

Proper dimensioning of the premixer and combustor place a restriction on the transverse (x) 
cell dimension in the computational mesh. This restriction forces the cross-sectional area 
ratio of the mesh to be slightly higher than the corresponding experimental value. (Note: 
The KTVA-n computer program, in original form, does not accommodate a symmetry plane 
in Cartesian meshes. Appendix B shows the four places in the code where modifications are 
made so that a symmetry plane can be utilized in this analysis.) 

3.2 Boundary Conditions 
Solid Boundaries 

At the walls of the dump combustor, normal components of velocity are set to zero. To 
determine the tangential velocity components, the turbulent boundary layer is considered 
in two parts, separated by a critical Reynolds number R c . The first part is the laminar 
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sublayer adjacent to the wall whose non-dimensional velocity profile is defined as: 

4- = t 1/2 (?s* e ) • 

u 

£ is the Reynolds number based on gas velocity u relative to the wall evaluated a distance 
y from the wall: 

y _ P yu 
p - (T) 

u*, called the wall-friction velocity, results from the dimensional analysis of 
u = u (p, x , p, y) where x w is the wall shear stress. Its value is determined by: 



In the second part of the two-part boundary layer model, the turbulent overlap and outer 
layers are combined. Here, the law-of-the-wall formulation is used in conjunction with the 
one-seventh power law to define the non-dimensional velocity profile: 

4 = ita(c,„? 7/8 )+B (S>jy 

u 

where c lw = 0.15, B = 5.5, and Karmann constant k = 0.4327. The critical Reynolds 
number R c , determined by equating the two expressions above for is 114. 

Isothermal walls are chosen for the dump combustor model based on the experimental setup 
of Proctor and T’ien. The wall heat flux using the two part boundary layer treatment 
above, is: 

j 1/(^4 r) ^<R C ) 

p« c p (T-T w ) ^ pr 4r+ (C>^ c ) 

where Pr t is the fluid Prandtl number in the laminar sublayer and T w is the wall 
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temperature. 

Symmetry Plane 

At the symmetry plane, a free slip condition is imposed on tangential components of 
velocity while normal components are set to zero. The symmetry plane is adiabatic and 
impervious to the transport of mass, momentum, and turbulent kinetic energy. 

Inflow - Outflow 

Beyond the inflow and outflow boundaries, fluid transport is assumed isentropic. Species 
densities and mixture pressure and temperature at the inlet and outlet planes are related to 
corresponding stagnation quantities using the equation of state and the isentropic relation: 

p_o _ rM 1 '' 1 ' 

P IP J 

where y is the ratio of specific heats. Turbulence quantities at the inflow and outflow 
boundaries are also related to corresponding stagnation values. 

33 Initial Conditions 

Species densities and mixture pressure and temperature are prescribed using the equation 
of state in three separate regions: 1) upstream of the premixer inlet; 2) the dump combustor 
flow field; and 3) downstream of the combustor exit. Dry-air (0 2 +3.76N 2 ) exists at rest in 
the combustor at the beginning of each program run. 
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CHAPTER 4 

COMPUTATIONAL PROCEDURE 

4.1 Start-up 
Cold Flow 

A steady state cold flow calculation is made first (Figure 6). This is achieved quickly since 
the geometry of the mesh is simple. An important note regarding this and further 
calculations is that the pressure gradient scaling (PGS) method provided with the KIVA-II 
code is not used. The purpose of this method is to improve the efficiency in the calculation 
of low speed flows by narrowing the difference between the acoustic and convective time 
scales so that the Courant sound speed stability condition is more easily satisfied. This is 
achieved by lowering the sound speed in the fluid by a factor cc, called the PGS parameter. 
Though the PGS method improves the efficiency of low speed implicit calculations 
significantly, it sacrifices accuracy in the modeling of acoustic disturbances. Since such 
disturbances are an important aspect of the combustion instability mechanism described in 
section 1.3, the PGS method is not implemented. 

The Reynolds number based on the hydraulic diameter of the premixer varies with the 
velocity, temperature, and equivalence ratio of the incoming mixture. Using the kinematic 
viscosity of air at the inlet temperature, the Reynolds number is on the order of l(r for all 
inputs of velocity and equivalence ratio in this study, which corresponds to a turbulent flow 

regime. 

Fuel Injection and Ignition 

The existence of the exit plate at the end of the mesh makes ignition difficult because it 
creates a semi-confined flow field that often induces large amplitude pressure fluctuations 



during ignition which, in turn, inhibit the formation of a stable flame on the comer of the 
step. In the experiment this problem is avoided by turning on the igniter at the same time as 
the fuel is injected and then increasing the equivalence ratio to a stable flame regime in no 
less than 30 seconds [Proctor, personal correspondence]. This prevents a buildup of fuel in 
the combustion chamber prior to and during ignition. Unfortunately, performing a 
calculation to 30 seconds requires excessive computing time even on a Cray-XMP. To 
circumvent this difficulty, the pre-exponential factor in the Arrhenius expression for the fuel 
oxidation reaction is linearly increased in time. This control of the energy release from 
combustion is a means of avoiding rapid temperature increases that, in turn, cause extreme 
pressure changes. 

Fuel injection is accomplished by changing the stagnation densities upstream of the inlet 
plane so that propane is added to the dry-air mixture at a specified equivalence ratio. A cold 
flow steady state velocity field comprising only dry-air is always achieved prior to fuel 
injection and ignition. (Note: The quick equilibrium solver CHMQGM in KIVA-II requires 
trace amounts of H 2 0 and C0 2 in every cell to avoid divide-by-zero errors. These 
compounds are included at partial densities on the order of 10' 25 in the cold flow mixture 
to fulfill this requirement.) 

Ignition is approximated in the computational model by adding specific internal energy to 
the fluid in a preselected group of cells for a specified time interval. The rate of heat addition 
is controlled by an input parameter. It is found that a four percent increase per time step 
(compounded every time step) produces a successful ignition for the equivalence ratios 
considered in this study (0.4 < 0 < 0.7). 

A fourth-order node coupler is provided with the KIVA-II code to smooth velocity 
fluctuations that may arise during the Lagrangian portion of the ICEd-ALE solution 
procedure. Though the node coupler is not needed to achieve a good steady state cold flow 
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solution, it is used during the combustion phase of the calculation to smooth velocity 
fluctuations that appear near the exit. 

4.2 Calibration of the Chemical Kinetics Model 

Proctor and Tien find that a stable flame attached to the dump plane at the end of the 
premixer occurs under the following conditions: 

(1) inlet air temperature - 850K 

(2) premixer wall temperature - 750K 

(3) average premixer velocity - 2133 cm/sec 

(4) equivalence ratio - 0.44 

(5) Reynolds number based on premixer channel width - 6562 1 

As the equivalence ratio is raised from 0.44, they observe an increase in the amplitude of 
combustor pressure oscillations. A maximum amplitude is reached just prior to flashback 
at an equivalence ratio of 0.56 (Figure 5a+b). 

In order to achieve flashback under the same conditions in this numerical analysis, it is 
necessary to calibrate the chemical kinetics model by varying the pre-exponential factor in 
the Arrhenius expression for the reaction rate coefficient of the fuel oxidation reaction. For 
the combustion of C 3 H 8 , Westbrook and Dryer [15] specify a pre-exponential factor of 
8.6X10 11 , noting that this value yields flame speeds that correlate well with experimental 
data. They also note, however, that global reaction rate mechanisms can vary with the type 
of experiment or combustion environment. In this analysis, a stable flame cannot be 
achieved under the conditions cited above with a pre-exponential factor of 8.6X10 1 1 . 2 Using 
smaller values with the same order of magnitude, a stable flame is established. Then, upon 

1 . Using kinematic viscosity of air at 800K, 1 01 .3 kPa. 

2 . It is noted in Chapter 5 that when using the quasi-second order upwind differencing scheme provided 
with the KIVA-II code, a stable flame can be calculated with the pre-exponential factor equalling 
8.6x10’ ’. However, for reasons also described in Chapter 5, the partial donor-cell differencing scheme 
is used in nearly all computations. 
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raising the equivalence ratio, it is found that the occurrence of flashback is strongly 
dependent on the pre-exponential factor. This is to be expected since the fuel oxidation 
reaction is the dominate source of heat in the kinetics model and that the pre-exponential 
factor Strongly influences the rate of heat release. A value of 1.9xlO n is chosen as a 
benchmark pre-exponential factor since it yields a stable flame at an equivalence ratio of 
0.44 (Figure 7) and also elicits flashback at approximately 0.57, in good agreement with the 
Proctor and T’ien experiment. 1 



As mentioned in the previous section, the pre-exponential factor is linearly increased in 
time during ignition. Prior to this the pre-exponential factor was held fixed and this caused 
a severe pressure spike in the combustor that immediately drove the flame into the premixer. 
Figure 8 shows a more moderate pressure rise brought on by the rapid increase in the rate 
of combustion as the bulk fluid temperature rises. Note that the stable flame calculation 
exhibits acoustic oscillations which are most likely a result of the ignition. (That these 
oscillations are acoustic is shown in Appendix C). 

Premature flashback is also induced if the equivalence ratio is raised too quickly from the 
stable flame calculation of Figure 7. In the interest of computational efficiency, the time rate 
of increase in equivalence ratio is chosen to be 0.1 per second, a value slightly lower than 
the minimum necessary to avoid influencing the flashback result. Any further decrease in 
the rate of equivalence ratio ramping does not affect the occurrence of flashback. 


4.3 Simulations 


Beginning from the stable flame calculation of Figure 7, two independent computations are 


1. At this point in the analysis, a run is made to determine if the stable flame calculation of Figure 7 is 
mesh dependent. The aspect ratio (lengthAvidth) of the cells is changed from 2:1 to 1:1 making the 
mesh more coarse, and there is no difference in the result over that appearing in Figure 7. 

2. A calculation is made using a value of 0.05 per second for the rate of increase in the equivalence ratio, 
but there is no change in the flashback result over that when using a value of 0.1 per second. 
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performed. First, the equivalence ratio is linearly increased and second, the inflow velocity 
is linearly decreased. Since the inlet air temperature and premixer wall temperature are 
found by experiment not to be important parameters in the development and sustenance of 
flashback, they are not varied in this study. The goal is to: 

(1) investigate pressure oscillations and other aerothermochemical conditions in 
the flow field prior to flashback; and, 

(2) investigate vorticity as the flame becomes unstable. 

The flashback observed by Proctor and T’ien as a result of raising the equivalence ratio has 
been discussed in the previous section with reference to Figures 5a+b. Proctor and T’ien 
also induce flashback by decreasing the average premixer velocity from 2621 cm/sec to 
23 1 6 cm/sec when the inlet air temperature is 830K, the premixer wall temperature is 750K, 
and the equivalence ratio is held fixed at 0.48. In the interest of computational efficiency, 
the linear decrease of the inflow velocity in the numerical calculation begins from the stable 
flame result of Figure 7 rather than computing a new stable flame using the experimental 
conditions cited above. This approximation is deemed valid, particularly since the inlet air 
temperature and the premixer wall temperature are observed not to significantly influence 
the occurrence of flashback. 

Calculated results occur in a variety of useful forms: a) contour plots of pressure, 
temperature, and equivalence ratio in the flow field at specified times, b) time history plots 
of pressure, temperature, and equivalence ratio at prescribed locations in the flow field, c) 
streamline plots of the flow field downstream of the step and, d) printed output of 
aerothermochemical data for all computational cells. 
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CHAPTERS 

DISCUSSION OF RESULTS 
5.1 Increased Equivalence Ratio 

Figure 9 shows a calculated flashback sequence which is generated by linearly increasing 
the equivalence ratio from the stable flame value of O = 0.44. A characteristic feature of 
impending flashback is the partitioning of the hot, diffuse region downstream of the high 
temperature gradient ‘cone’ (Figure 9a). This is followed by a collapse of the ‘cone’ into 
the premixer (Figure 9b-i). Note that this sequence occurs very rapidly (approx. 0.04 
seconds), as expected. Also shown in Figure 9 is a sequence of nine streamline plots 
representing the same elapsed times in the calculation as the nine temperature contour plots. 
A disappointing feature of the streamline plots is the absence of coherent vortices in the 
turbulent shear layer extending from the dump plane. The occurrence of these vortices is 
predicted in step (3) of the combustion instability mechanism described in section 1.3. Of 
course, without this step, unsteady heat release oscillations do not develop and the self- 
excited, ‘closed loop’ character of the mechanism can never be sustained. 

Figure 10 shows the combustor pressure history during flashback which includes the 
sequence of Figure 9. The amplitude of the acoustic disturbances increases approximately 
30% 1 as the equivalence ratio is raised from 0.44 and then becomes irregular just prior to 
flashback. As flashback occurs a pressure jump is seen in the combustor, as expected. The 
amplitude of the pressure oscillations in this calculated result is substantially lower than 
during the onset of flashback in the Proctor and T’ien experiment (Figure 5b). It is possible 
that the higher amplitude oscillations seen in the experiment are a manifestation of the 
growth and convection of large-scale coherent vortices in the shear layer. Combustion 


1. This percentage is calculated by comparing the change in amplitude of the acoustic oscillations from 
the stable flame regime of Figure 8 to that of impending flashback observed in Figure 10. 
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within the vortices may then have caused unsteady heat release oscillations as described in 
step (4) of the instability mechanism in section 1.3. This could have been the driving 
mechanism for the flame flickering mode that Proctor and T’ien observe prior to flashback. 
Note that Proctor and T’ien do not propose this as a driving mechanism for the combustion 
instabilities they observe since their work predates that of Smith and Zukoski [9], Schadow, 
et.al. [11], Logan, et.al. [12], and others [18-20] which has been essential to the evolution 
of the instability mechanism of section 1.3. 

As a possible explanation for the aforementioned disparity in the calculated and 
experimental pressure oscillations prior to flashback, it should be noted that Proctor and 
T’ien raise the equivalence ratio from 0.44 up to the flashback limit of 0.60 in an elapsed 
time of a few minutes [Proctor, personal correspondence] while the calculated result has d> 
raised from 0.44 to 0.57 in approximately 1.1 seconds. This being so, there may not be the 
opportunity for coherent vortices to form in the shear layer of the computational model. As 
mentioned at the end of section 4.2 the rate of equivalence ratio ramping is lowered from 
0.1 per second (the value used to calculate the flashback sequence of Figure 9) to 0.05 per 
second, but there is no change in the solution. To determine whether enough time elapses 
in the numerical simulation for vortices to develop, it is instructive to examine the time 
scale for vortex formation in the region of the shear layer just downstream of the dump 
plane. First, consider a typical vortex having a 1 .0 cm diameter in this initial part of the 
shear layer. In addition, assume an average velocity of 1500 cm/sec in the vicinity of this 
developing vortex. 1 The time scale for formation should be the circumference of the vortex 
divided by the surrounding fluid velocity which is approximately 0.002 seconds. Based on 
this rough calculation, there appears to be sufficient time in the numerical calculation for 
large-scale vortices to develop in the shear layer. 


1. Examination of the printed output from the stable flame calculation shows that at a point 
approximately 4.0 cm downstream of the dump plane, the shear layer width is about 1 .5 cm and the 
average velocity across the shear layer is about 1500 cm/sec. 
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To substantiate the accuracy of the flashback calculation described in this section, a run is 
made in which the convergence criteria on pressure and temperature are lowered a whole 
order of magnitude (1 x lO -4 to 1 x lCf 5 and 1 x 10 -3 to 1 x lO -4 , respectively). This 
change yields no variation in the solution over that appearing in Figure 9. 

So, this begs the question of what is causing the flashback in the calculated result. Is it a 
consequence of the approximately 30% increase in the amplitude of the acoustic 
oscillations? It should be noted here that the computed flashback result of Figure 9 is 
achieved using a partial donor-cell differencing scheme. When the quasi-second order 
upwind (QSOU) differencing scheme is used, the combustor and premixer pressure 
oscillations are completely damped out and no flashback occurs for values of O up to 
approximately 0.7, beyond which no calculation is performed. A second run is made using 
the QSOU differencing scheme with a pre -exponential factor of 8.6 x 10 1 1 , consistent with 
the Westbrook and Dryer one-equation global reaction rate mechanism for propane 
oxidation. Then the equivalence ratio is raised at a rate of 0.1 per second up to 
approximately 0.70 and, once again, neither pressure oscillations nor flashback occur. So, 
it is clear that the calculated flashback result of Figure 9 does depend on the differencing 
scheme used. For this computation, the QSOU differencing scheme is clearly more 
diffusive even though it is formally more accurate than the partial donor-cell method [3]. 
An implicit and otherwise unsubstantiated conclusion can be made from this dependence 
on the differencing scheme: since flashback never occurs (up to O = 0.70) when acoustic 
pressure oscillations are damped out and does occur when there is a small increase in the 
amplitude of the oscillations, then perhaps this slight intensification of the acoustic field is 
the driving mechanism for flashback in the computed solution. 

Regardless of whether this conclusion is valid or not, there does not appear in the calculated 
result the large, coherent vortices and the flame flickering mode observed in experiment 
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prior to flashback. Consequently, the complete instability mechanism described in section 
1.3 is not accurately modeled in this calculation. As is mentioned in the next section in 
conjunction with decreasing the inflow velocity, this deficiency may be due in part to the 
spatial resolution of the premixer boundary layer and/or the k - e turbulence model. 

5.2 Decreased Inflow Velocity 

As outlined in section 4.3, a second part of this analysis is to linearly decrease the inflow 
velocity while holding the equivalence ratio fixed at 0.44. Two separate calculations are 
made; one using the partial donor-cell differencing scheme and the other using the QSOU 
differencing scheme. It is found that, as the inflow velocity decreases, the high temperature 
gradient ‘cone’ seen in Figure 7 shortens and the combustor pressure decreases linearly, but 
no flame flickering or flashback occur. In each calculation, the inflow velocity is decreased 
at a rate of 750 cm/sec 2 to 1050 cm/sec, which is approximately half the original inflow 
velocity of 2133 cm/sec. Figure 11 presents the final segment of the combustor pressure 
history for the calculation using partial donor-cell differencing. It shows that there is no 
increase in the amplitude of the acoustic oscillations over that appearing in the latter part of 
Figure 8. In addition, it is interesting to note that the oscillations have become more 
irregular, an effect most likely caused by the increased influence of the step on the 
propagation of the acoustic waves. 

Experimentally, when the inflow velocity is varied over a broad range, many turbulent 
shear layer instability frequencies result. This increases the likelihood that a constructive 
phase relationship will arise between the acoustic field and the shear layer instabilities as 
described in steps (l)-(3) of the instability mechanism in section 1.3. In an attempt to 
calculate a coupling between acoustic and shear layer oscillations, an additional 
computation is performed in which the inflow velocity is linearly increased from the stable 
flame result of Figure 7. Once again, two separate runs are made, one using partial donor 


cell differencing and the other .using QSOU differencing. In each calculation, the inflow 
velocity is increased at a constant rate of 1500 cm/sec 2 to a maximum of 4300 cm/sec, 
which is approximately twice the original inflow velocity of 2133 cm/sec. It is found that 
the high temperature gradient ‘cone’ elongates, as expected, and the combustor pressure 
increases linearly, but again no flame flickering or flashback occur. 

Following are two possible explanations for the absence of combustion instabilities in the 
calculated result as the inflow velocity is varied. First, as is noted in step (2) of the 
instability mechanism in section 1.3, the velocity profile and flow structure of the premixer 
boundary layer affect the occurrence of combustion instabilities. The turbulent boundary 
layer in the premixer is thin and is therefore not spatially resolved by the boundary cells in 
the premixing region of the computational mesh. The law-of-the-wall boundary layer 
model describes the effects of turbulence on the wall shear stress and heat transfer 
coefficient but does not describe the turbulent boundary layer vortices themselves. The 
development and transport of these vortices may provide the seed for the development of 
coherent vortices in the shear layer downstream. To resolve turbulent boundary layer 
vortices, the mesh could be refined in the vicinity of all solid boundaries, but the expense 
of such a calculation would be extreme. 

Before undertaking such an expensive calculation, a second possible explanation for the 
absence of vortices in the shear layer of the calculated result should be considered. The 
k-e model has existed for some time [14] and it often does quite well to calculate a local 
value for the effective turbulent viscosity based on solutions to transport equations for 
turbulent kinetic energy and its dissipation rate. However, it has already been shown that 
the effects of numerical diffusion can play an important role in the occurrence of flashback, 
which indicates that this is a diffusion-sensitive phenomenon. The magnitude of turbulent 
diffusion may play a similarly important role in permitting/suppressing development of 
coherent vortices from shear layer instability. Development of a new turbulence model or 
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even refinement of the Jc — £ .model parameters are beyond the scope of this work. 
Turbulence modeling is a difficult field and, despite the efforts of some of the best applied 
mathematicians, gains are made only very slowly. 






. CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 
6.1 Conclusions 

A computed flashback sequence is achieved using the KIVA-II code by linearly increasing 
the equivalence ratio from the stable flame calculation of Figure 7. It is observed that the 
following two synergistic features that are identified experimentally by Proctor and T’ien 
and others do not appear in the calculated result: 

(1) the conception and growth of coherent vortices in the turbulent shear layer, 
and; 

(2) high amplitude combustor pressure oscillations coincident with flame 
flickering prior to flashback. 

The absence of these features elicit the conclusion that the computational model does not 
accurately reproduce the experimental instability mechanism described in section 1.3. It is 
noted that the near- wall grid resolution and/or the k-z turbulence model may be 
inadequate to simulate step (2) of this mechanism. This would then prevent the coupling of 
acoustic and shear layer frequencies as required for the conception of coherent vortices and, 
in tum, high amplitude combustor pressure oscillations. Also noted is the disparity in the 
time elapsed between the Proctor and T’ien experiment and the numerical calculation as the 
equivalence ratio is raised from 0.44 to the flashback limit. The question is raised as to 
whether vortices can form during the much shorter elapsed time of the numerical 
simulation. It is shown by means of a rough calculation of the time scale of vortex formation 
in the shear layer just downstream of the dump plane that sufficient time exists for the 
development of coherent vortices. 

As mentioned in Chapter 5, the computed flashback result does depend on the finite 
differencing scheme used. The QSOU differencing scheme is more diffusive in this 
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calculation since it damps out {he acoustic oscillations that occur when the partial donor- 
cell method is employed. That flashback occurs only when the acoustic oscillations are 
present leads to the conclusion that the intensification of the acoustic field resulting from 
the linear increase of the equivalence ratio may be the driving mechanism for the calculated 
flashback. 

The inlet velocity is linearly decreased from the stable flame calculation of Figure 7 using 
both partial donor-cell and QSOU differencing, but neither an increase in the amplitude of 
the acoustic pressure oscillations nor flashback occur. In an attempt to induce coupling 
between the acoustic and shear layer oscillations as described in steps (l)-(3) of the 
instability mechanism in section 1 .3, the inflow velocity is also increased over a broad range 
but no excitations are calculated. These results are attributed to the course spatial resolution 
of the premixer boundary layer and to uncertainty in the ability of the k - e model to 
adequately simulate the small-scale turbulence in this boundary layer. As suggested 
previously, these features of the numerical model may preclude the development of 
turbulent shear layer instability waves which, in turn, would prevent the coupling of 
acoustic and shear layer oscillations necessary to form coherent vortices. 

6.2 Recommendations for Further Analysis 
Exit Aperture and the Acoustic Field 

As noted in section 3.1, the exit aperture dimensions in the computational mesh are chosen 
so that the exit/combustion chamber cross-sectional area ratio correlates with the 
experiment. This properly models the incompressible flow resistance, but what of the 
acoustic field in the entire dump combustor? 

The standing wave patterns comprising an acoustic field develop from the interference of 
longitudinally traveling waves that are reflected at each end of the dump combustor. The 
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nature of the reflections depends on the open-to-closed area ratios of the inlet and exit. As 
an illustration, consider a dump combustor in which the inlet geometry is held fixed while 
the exit dimensions are varied. When the end of the combustion chamber is completely 
closed (exhausting at another location), the fluid elements adjacent to the end wall are 
stationary (no slip at the wall), giving rise to an acoustic node. Conversely, the 
establishment of a completely open boundary at the end of the combustor would allow for 
nearly unrestricted movement of the fluid elements along the exit plane and an approximate 
antinode would occur. (A descriptive analogue is to consider the difference in the wave 
patterns for two longitudinally vibrating beams, one with both ends fixed and the other 
having one end free). Between these two extremes of fully closed to fully open are infinitely 
many exit geometries, each associated with a unique acoustic field. Therefore, the 
frequency of the standing wave fundamental and harmonic modes are a function not only 
of the length of the dump combustor, but also the inlet and exit dimensions. 

In this analysis, the inlet geometry is easily specified since the experimental premixer is 
purely two-dimensional. Conversely, the geometry of the two-inch hole exhausting the 
experimental combustion chamber is three-dimensional and therefore more difficult to 
model acoustically. An even greater task, if at all possible in a two-dimensional analysis, is 
to select an exit geometry that simultaneously elicits the correct flow resistance and 
acoustic response. A compromise is made in this analysis by ensuring proper flow 
resistance while risking inaccuracy in modeling the acoustic response of the exit. Since the 
acoustic field is an important mechanism in the occurrence of flashback, this compromise 
may cause the computational model to depart significantly from the experiment. 
Consequently, if the Proctor and T’ien experiment is to be used as a basis for further 
numerical analysis, it is recommended that a three-dimensional model be developed. A 
simpler tactic would be to choose an experiment that uses an apparatus which has two- 
dimensional features throughout such as that in reference [12] (Karagozian, personal 
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correspondence). Perhaps then .the CPU time and data storage that would have otherwise 
been consumed by computing in a three-dimensional volume could be utilized to refine the 
mesh in the boundary layer regions. 

Symmetry Plane and the Acoustic Field 

With the symmetry plane in the computational mesh, accurate simulation of transverse 
acoustic wave propagation and reflection is not possible in the numerical model. These 
acoustic waves are most likely present in the Proctor and T’ien experiment and may excite 
high frequency ‘screech’ modes of combustion instability. However, it is the lower 
frequency modes of instability induced by longitudinal acoustic wave propagation that have 
been the focus of recent experimental investigations [9,11,12,18-20] and which are 
described in the combustion instability mechanism in section 1.3. This being so, it is 
deemed worthwhile to continue to utilize a symmetry plane in further analysis. (Note: The 
asymmetries observed in Figure 5a are most likely a manifestation of asymmetric features 
of the Proctor and T’ien experimental rig such as the location of pressure transducers and 
stream thermocouples.) 

Streaklines vs. Streamlines 

Streamlines (loci of constant velocity gradient) are used to examine features of the velocity 
field, particularly the turbulent shear layer. In the review of this thesis, it was noted that the 
streamline plots may not reveal the large-scale coherent vortices associated with the 
combustion instability modes observed in experiment. Though not presented, velocity 
vector field plots are consistently used as an analytical tool to substantiate the streamline 
plots. In further analysis, however, it may be worthwhile to use plots of marker particle 
paths, or streaklines, since they may be more revealing than streamlines, particularly with 
regard to shear layer vortices. 
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Chemical Kinetics - One vs. Multi-Equation Model 

The one-equation global reaction rate mechanism for propane oxidation developed by 
Westbrook and Dryer [15] appears to work well in this numerical analysis. A more rigorous 
test of this mechanism would be to see how well it models combustion within developing 
coherent vortices in the turbulent shear layer as predicted in step (3) of the instability 
mechanism in section 1.3. Unfortunately, as described in the previous chapter, coherent 
vortices do not appear in the shear layer of the calculated result. To say that this is a 
manifestation of the approximate nature of the global reaction rate mechanism is probably 
an overstatement. In fact, as suggested in the previous section, the disparity between the 
features of the calculated flashback sequence and those of the instability mechanism in 
section 1.3 is most likely due to the inability of the computational model to adequately 
simulate the first three steps in the mechanism that describe the coupling of acoustic and 
turbulent shear layer frequencies. This phenomenon is, of course, convective in nature and 
therefore not solely dependent on the reaction rate. However, since the reaction rate 
mechanism controls the rate of heat release from the oxidation of the fuel and this, in turn, 
governs the flame shape and structure, it may be worthwhile upon further analysis to 
incorporate a more elaborate reaction rate mechanism. The analyst will need to decide how 
much accuracy is gained for the added computing time required by a multi-equation 
reaction rate mechanism. 
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Figure la: Experimental Dump Combustor after Logan, et.al. [12] 
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Figure lb: Map of Combustor Modes after Logan, et. al. [12] 
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Figure 2a: KIVA-II - Turbulent Flow over a Step 
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Figure 3a: Cinematographic Schlieren Record of Stable Combustion at O - 0.57, 

after Keller, et.al. 
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Figure 3b: KIVA-II Numerical Model of Stable Flame in Fig. 3a 










Figure 5a: Proctor and T’ien Experimental Flashback Sequence 



Different stages of burning, (a) and (b), lean 
flames, equivalence ratio, 0.41; (c) recirculation 
zones formed behind step, equivalence ratio, 0.44; 

(d) flame began flickering, equivalence ratio, 0.56; 

(e) maintained flashback, equivalence ratio, 0.60. 

Inlet air temperature, 850K; premixer wall temperature, 
750K; average premixer velocity, 70 ft/s; fuel injector 
8 in. upstream of step. 
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Figure 5b: Typical Combustor Pressure History 
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L - low contour line 

H - high contour line 
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Figure 6: Cold Flow Steady State Calculation 



Figure 7: Stable Flame Calculation at <l> = 0.44 







Figure 8: Combustor Pressure History of Reactive Flow Start-Up 
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Figure 11: Combustor Pressure History as Inflow Velocity is Decreased 
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> Appendix A 

Appendix A shows the procedure used to compute inlet stagnation species densities of the 
reactant mixture. Inlet stagnation pressure and temperature remain constant in this study so 
the only variable in the following procedure is equivalence ratio. Note that superscript o 
identifies stagnation properties. 

A stoichiometric reaction is first determined using a simple atomic balance: 

\ 

C 3 // 8 + 5 ( 0 2 + 3.76jV 2 ) 3C0 2 + 4// 2 C>+ 18.8N 2 . 

The equivalence ratio O is defined on a mass basis as: 

f m fuel ] 


<D = 


\ m oxidant ) ac t ua i 


m 


'fuel 


< m oxidant J st0lc h 


'.'-i 

¥4 




L 


Note that the oxidant is taken to be dry air (0 2 +3.76N 2 ). Using the stoichiometric reaction 
above and the molecular weights of C3H8, 0 2 , and N 2 : 


m /« e / ^ 
K. m oxidant Jstoich 


= 0.0641025 . 


The following expression is derived for the stagnation species density of 0 2 after defining 
a volume of mixture, V,^, at the inlet and applying the equation of state there: 


P n. = 


mix 


2 R T° ■ 

u mix 


. ,,,' V M<t>(0.0641025) _ 4.76 

1 + 3. /O — — I tt; + 


Wr 


W. 


c 3 h. 


W, 


-1 


The species densities of C 3 H 8 and N 2 are then determined using the following relations: 

Pc,//, = (P^+py^ (0-0641025) 
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p k ■ 3 - 76 «J 


The expressions above for p^, Pc 3 // 8 ’ 316 incorporated into a'subroutine called 

FUELAIR (shown below) which is added to the KIVA-II code for this analysis. Note that 
in its original form KIVA-II uses variables ‘pamb’ and ‘tempamb’ to identify stagnation 
states upstream of the inlet and downstream of the exit. For this analysis, variables ‘pinO’ 
and ‘tinO’ are substituted for ‘pamb’ and ‘tempamb’ where the inlet stagnation state is used 
in computation. 


*deck fuelair 

subroutine fuelair 
c 


c 

c 

c 

c 

c 

c 

c 

c 

*call 
c 
c 
c 

c === 
c 


subroutine fuelair computes the stagnation densities 
of a c3h8-o2-n2 mixture using input equivalence ratio, 
stagnation pressure, and temperature. 

fuelair is called by: chem 


comd 


the variable 'stoich' is the stoichiometric ratio of 
fuel to oxidant as defined by reaction below: 

c3h8 + 5 (o2 + 3 . 76n2 ) > 3co2 + 4h2o + 18.8n2 

note: oxidant is assumed to be dry air (o2 + 3.76n2). 
stoich= ( 1 . *mw( 1 ) )/( (5.*mw(2) )+(18 . 8*mw( 3 ) ) ) 


alpha=( 1 . +3 . 76* (mw( 3 )/mw( 2 ) ) ) 

& *( (eqrin*stoich)/mw(l) )+( 4 . 76/mw( 2 ) ) 

spdinO ( 2 ) = ( pinO/ ( rgas*tinO ) ) /alpha 
spdinO ( 3 ) =3 . 76*spdin0( 2) *(mw( 3 ) /mw(2) ) 
spdinO ( 1 ) = ( spdinO ( 2 ) +spdinO ( 3 ) ) *eqrin*stoich 
c 

c === recompute roinO, rgamin, and dprefin which were 
c === originally computed in subr. setup, 
c 


it=int( 0 . 01*tin0) 

elo=0 . 0 

ehi=0 . 0 

roin0=0 . 0 

do 5270 isp=l,nsp 

elo=elo+spdinO ( isp) *ek ( it+1 , isp) 

ehi=ehi+spdinO ( isp) *ek ( it+2 , isp) 

roin0=roin0+spdin0 ( isp) 

5270 continue 

csubv=0 . 01* ( ehi-elo) /roinO 

rgamin=l .0/(1. 0+pin0/(roin0*csubv*tin0) ) 

dpref in=botin*0 . 5*roin0*win**2 

return 

end 
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Appendix B 

Appendix B shows segments of four subroutines in KIVA-II that required modification so 
that a symmetry plane could be used in the computational model. A simple change in the 
logical argument of each ‘if’ statement below treats the symmetry plane as if it were a 
symmetry axis in polar coordinates: adiabatic and impervious to mass, momentum, and 
turbulent kinetic energy transport. Note that the symmetry plane is located along the first 
radial plane (i=l). 

Subroutine BCEPS: 

c =20 if ( ( i . eq . 1 . and. cyl . eq . 0 . ) .or. (i.gt. 1 .and.f (i4-l) . eq.O . ) ) then 
20 if ( ( i . eq. 1 .and. cyl .eq. 1 . ) .or. ( i .gt . 1 . and. f ( i4-l ) . eq . 0 . ) ) then 
if 4=i4 
go to 30 

elseif ( (i . eq.nx .and. rtout . eq. 0. ) .or . ( i . It .nx .and. f (il) .eq.O. ) ) 


Subroutine BCRESE: 


c 


ef acik=cvmgt (0 .0,efacik,i.eq.l 
ef acik=cvmgt ( 0 . 0 , ef acik, i . eq. 1 
ef acik=cvmgt ( 0 . 0 , ef acik , i . gt . 1 
ef acik=cvmgt ( 0 . 0 ,ef acik, i . eq.nx 
ef acik=cvmgt ( 0 . 0 ,ef acik, i . It . nx 


•and. cyl.eq.0.0) 
•and. cyl.eq.1.0) 
•and. f (i4-l) .eq.0.0) 
. and . rtout .eq.0.0) 
•and. f(i4+l) .eq.0.0) 


Subroutine LAWALL: 

c =40 if ( (i.eq.l. and. cyl. eq.O. ) .or. ( i -gt . 1 . and. f ( im jk) . eq . 0 . ) ) goto 50 
40 if ( (i.eq.l. and. cyl. eq.l. ) .or. (i.gt.l.and.f (imjk) .eq.O. ) ) goto 50 
if ((i. eq.nx .and. rtout . eq. 0 . 0) .or. 

& (i.lt.nx .and. f (il ) . eg. 0 . 0) ) go to 60 
go to 80 


Subroutine SETUP: 

c === if ((i.eq.l .and. cyl. eq.O.) .or. 
if ((i.eq.l .and. cyl. eq.l.) .or. 

& (i.gt.l .and. f ( i4-l ) . eq. 0 . ) ) go to 5110 

if ( icart3 . eq. 0) go to 5140 

if(j.eq.l .or. f ( i4-nxp) . eq. 0 . ) go to 5120 
if(j.eq.ny .or. f ( i4+nxp) . eq. 0 . ) go to 5130 
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Appendix C 


The stable flame calculation of Figure 7 exhibits acoustic oscillations that are seen in the 
latter part of Figure 8. As noted in section 5.1, these oscillations increase in amplitude as 
the equivalence ratio is raised linearly in time. Figure 10 shows the acoustic disturbances 
as they become irregular prior to flashback. To show that these oscillations are acoustically 
driven, the following calculation is made: 

Assume the fluid in the dump combustor is air at 1200K. The speed of sound c in the fluid 
is then: 

r 'tia't * io^ 

c = JyRT = ( 1.4) ( 28 97 -- "- ) ( 120 °) =70, OOOcm/sec 


where the value 28.97 is the molecular weight of air which is needed for proper unit 
cancellation. 


Define Fas the period of the regular, symmetric oscillations seen in Figures 7 and 10: 


0.05 sec 
19cycles 


0.0026 sec . 


Note that an accurate reading of the number of cycles occurring in 0.05 seconds is obtained 
from pressure history plots representing the time interval between Figures 7 and 10. In the 
interest of brevity, these plots are not presented. 


The overall length of the dump combustor is 43.18 cm. Define t* as the time elapsed as a 
sound wave travels this distance: 

* _ 43.18cm 

* 70, OOOcm/sec 

The inlet in the computational model is velocity specified and can therefore be considered 
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an acoustic node. Conversely, the outlet, being pressure specified, is considered an 
approximate antinode. Consequently, a quarter wave exists in the computational mesh so 
that t* must be multiplied by 4 to determine the period of acoustic oscillations, defined as 
x below: 

x = 4 x f* = 0.0025 sec . 


Finally, comparing T with x, it can be seen that the oscillations are acoustic. 
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Appendix D 


The input file used in the computations relating to the Proctor and T’ien experiment is listed 
below. There are a few variables that are added to the original KTVA-II input file for this 
analysis. Shown first is a description of these supplemental variables and then the input 
listing. 


nctap9 

nclst 

fturb 


tlfuel 

winrate 

eqrate 

pinO 

spdinl-12 


cycle interval between writes to post-processing file tape9. 
beginning cycle for history plot arrays. 

multiplicative factor used to determine the influence of effective turbulent 
viscosity on the reactive flow calculation. Appears only in subr. VISC: 
terme=fturb*cmu*tke( i4) **2*ro(i4 ) /epsden 
where terme is the effective turbulent viscosity term in the equation for 
amu(i4), the cell-specific bulk fluid viscosity, 
time at which fuel is fluxed in with dry air at the inlet plane, 
rate of linear increase in inflow velocity, win. 
rate of linear increase in inflow equivalence ratio, eqrin. 
stagnation pressure upstream of the inlet plane, 
species stagnation densities upstream of the inlet plane. 


Also, as noted in section 4. 1 , the pre-exponential factor cf(l) for the fuel oxidation reaction 

is linearly increased in time during ignition. This is accomplished at the beginning of 

subroutine CHEM as follows: 

if (t .gt .tlign .and. cf ( 1) . It. 1.9e+ll) 

& cf (l)=(1.9e+ll/0.60)*(t-tlign) 

Activation temperature ef(l), temperature exponent zetaffl), and species concentration 
exponents ae(isp,l) are obtained from reference [15]. 
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maine k053190 16x1x69 2-d proctor, w/ pdc, chemeq 053190 



irest 

25 


ipost 

1 

? ; 

nx 

16 


ny 

1 


nz 

69 


lwall 

1 


nchop 

0 


lpr 

1 


jsectr 

0 


irez 

0 

— 

ncf ilm 

99999 


nctap8 

99999 


nctap9 

99999 



nclst 

45317 


nclast 

99999 


caf ilm 

9 . 99e+9 


caf in 

9 . 99e+9 


cadump 

9 . 99e+9 


dcadmp 

9 . 99e+9 

’•Tj i 

angraom 

0 . 0 

— 

cyl 

0.0 


dy 

0.5 

’1* tn:.: 

pgssw 

0.0 


sampl 

0.0 

r 

dti 1 

. . 04167e-4 

- - 

dtmxca 

9 . 99e+9 

v - e i 

dtmax 

9 . 99e+9 


tlimd 

3600.0 

- 

twf ilm 

9 . 99e+9 

M 

twf in 

9 . 99e+9 

_:7: 

fchsp 

0.25 

v&i- 

stroke 

9.55 

* 

"'14- F 

squish 

0.181934 


rpm 

0.0 

SI 

atdc 

0.0 


conrod 

16.269 

- > 

offset 

0.0 


swirl 

0.0 


swipro 

0.00 

ly f 

thsect 

360.0 

x^-'st r 

epsy 

1.0e-3 


epsv 

1.0e-3 

:"7 n : 

epsp 

1 . Oe-4 

fe j 

epst 

1.0e-3 


epsk 

1 . Oe-3 


epse 

1.0e-3 

* p 

gx 

0.0 


gy 

0.0 


gz 

0.0 


tcylwl 

750.0 


thead 

750.0 


tpistn 

750.0 


tvalve 

: 750.0 



tempi 

800.0 

pardon 

1.0 

aO 

0.0 

bO 

1.0 

anc4 

0.05 

adia 

0.0 

anuO 

0 . 0 

fturb 

1.0 

visrat- 

. 66666667 

tcut 

900.0 

tcute 

1200.0 

epschm 

0.02 

omgchm 

1.0 

tkei 

4 . 55e+4 

tkesw 

1 . 0 

sgsl 

0.0 

uniscal 

0.0 

airmul 

1 . 457e-5 

airmu2 

110.0 

airlal 

252.0 

airla2 

200.0 

expdif 

0.6 

prl 

0.74 

rpr 

1.11 

rprq 

1.0 

rpre 

0.769231 

rsc 

1.11 

xignit 

1 . 50e+3 

tlfuel 

2 . 20e-l 

tlign 

2 . 20e-l 

tdign 

0.20 

calign 

9 . 99e+9 

cadign 

9 . 99e+9 

iignll 

1 

iignrl 

2 

jignf 1 

1 

jigndl 

1 

kignbl 

32 

kigntl 

32 

iignl2 

0 

iignr2 

0 

jignf 2 

0 

jignd2 

0 

kignb2 

0 

kignt2 

0 

kwikeq 

1 

numnoz 

1 

numvel 

1 

in jdist 

1 

kolide 

0 

tlinj 

-1.0 

tdin j 

-1.0 

calin j 

9 . 99e+9 

cadin j 

9 . 99e+9 
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tspmas 

0.0116 

pulse 

1.0 

tnparc 

0.0 

rhop 

0.7436 

tpi 

350.0 

turb 

0.0 

breakup 

0.0 

evapp 

0.0 

drnoz 

0.0 

dznoz 

0 . 0 

dthnoz 

0.0 

tiltxy 

0 . 0 

tiltxz 

0.0 

cone 

62.5 

dcone 

12.5 

anoz 

1.0 

smr 

5 . 00e-4 

ampO 

0.0 


4000.0 


npo 17 

nunif 0 


1 

1 

0 . 00000 

0.0 

2 

1 

0.31750 

0.0 

3 

1 

0.63500 

0.0 

4 

1 

0.95250 

0.0 

5 

1 

1.27000 

0.0 

6 

1 

1.58750 

0.0 

7 

1 

1.90500 

0.0 

8 

1 

2.22250 

0.0 

9 

1 

2.54000 

0.0 

10 

1 

2.85750 

0.0 

11 

1 

3.17500 

0.0 

12 

1 

3.49250 

0.0 

13 

1 

3.81000 

0.0 

14 

1 

4.12750 

0.0 

15 

1 

4.44500 

0.0 

16 

1 

4.76250 

0.0 

17 

1 

5.08000 

0.0 

nho 


0 


square 

0.0 


rcornr 

0.0 


nstrt 

0 



icont 00001000000000000000000000 

mirror 1 

nvzone 0 

nvwec 0 

nvpvec 0 

nvcont 0 


nsp 

12 






c3h8 

rhol 

0.0 





o2 

rho2 

9 . 8544e-5 

mw2 

32.000 

htf 2 

0.0 

n2 

rho3 

3 . 2440e-4 

mw3 

28.016 

htf 3 

0.0 

co2 

rho4 

1 . 000e-25 

mw4 

44.011 

htf 4 

-93.965 

h2o 

rho5 

1 . 000e-25 

mw5 

18.016 

htf 5 

-57.103 


h 

rho6 

0.0 

mw6 

h2 

rho7 

0.0 

mw7 

o 

rho8 

0.0 

mw8 

n 

rho9 

0.0 

mw9 

oh 

rholO 

0.0 

mwlO 

CO 

* rholl 

0.0 

mwll 

no 

rhol2 

0 . 0 

mwl2 

rtout 

0 . 0 



topout 

1.0 



botin 

1.0 



distamb 

0.0 



pamb 

9 . 7500e+5 



tkeamb 

4 . 55e+4 



sclamb 

0.1905 



spdaml 

0.0 




spdam2 2.6455e-4 
spdam3 8.7086e-4 
spdam4 1.000e-25 


spdam5 

1 . 000e-25 

spdam6 

0.0 

spdain7 

0.0 

spdam8 

0.0 

spdam9 

0.0 

spdaml 0 

0.0 

spdaml 1 

0.0 

spdaml 2 

0.0 

win 

2133.0 

winrate 

1 . 5e+3 


1.008 

htf 6 

51.631 

2.016 

htf 7 

0.0 

16.000 

htf 8 

58.989 

14.008 

htf 9 

112.520 

17.008 

htf 10 

9.289 

28.011 

htf 11 

-27.200 

30.008 

htf 12 

21.456 


eqrate 0.025 
pinO 9.7500e+5 
spdinl 0 . 0 

spdin2 9.8544e-5 
spdin3 3.2440e-4 
spdin4 1.000e-25 


spdin5 l.OOOe- 
spdin6 0 . 0 

spdin7 0 . 0 

spdin8 0 . 0 

25 





spdin9 0 . 0 
spdinlO 0.0 
spdinll 0.0 
spdinl2 0.0 






nrk 4 






cfl 0.0 

efl 

1 . 5110e+4 zfl 

0.0 


cbl 0.0 

ebl 

0 

. 0 zbl 

0.0 


ami 1 5 

0 0 

0 

0 0 0 0 

0 

0 0 

bml 0 0 

0 3 

4 

0 0 0 0 

0 

0 0 

ael 0.100 1. 

650 0 

.000 

0.000 0.000 

0.000 

0.000 

0.000 0 . 

000 0 

.000 

0.000 



bel 0.000 0. 

000 0 

.000 

0.000 0.000 

0.000 

0.000 

0.000 0 . 

000 0 

.000 

0.000 



cf 2 1 . 5587el4 

ef 2 

6 . 7627e+4 zf2 

0.0 


cb2 7 . 5000el2 

: eb2 

0 

. 0 zb2 

0.0 


am2 0 1 

2 0 

0 

0 0 0 0 

0 

0 0 


0.000 

0.000 


68 


bin 2 

ooooo 

0 ' 

0 0 2 0 

ae2 

0.000 0.500 1.000 

0.000 0.000 0.000 


0.000 0.000 0.000 

0.000 

be2 

0.000 0.000 0.000 

0.000 0.000 0.000 


1.000 0.000 0.000 

1.000 

cf 3 

2.6484el0 ef3 5.9418e+4 

zf 3 1.0 

cb3 

1.6000e+9 eb3 1.9678e+4 

zb3 1 . 0 

am 3 

0 2 10 0 

0 

0 0 0 0 

bm3 

0 0 0 0 0 

0 

0 2 0 0 

ae3 

0.000 1.000 0.500 

0.000 0.000 0.000 


0.000 0.000 0.000 

0.000 

be 3 

0.000 0.000 0.000 

0.000 0.000 0.000 


0.000 0.000 0.000 

1 . 000 

cf 4 

2 . 1230el4 ef4 5.7020e+4 

zf 4 0.0 

cb4 

0 . 0 eb4 0 . 0 

zb4 0 . 0 

am4 

0 0 1 0 0 

0 

0 0 0 2 

bm4 

ooooo 

2 

0 0 0 0 

ae4 

0.000 0.000 0.500 

0.000 0.000 0.000 


0.000 1.000 0.000 

0.000 

be 4 

0.000 0.000 0.000 

0.000 0.000 1.000 


0.000 0.000 0.000 

1.000 

nre 

6 



asl 

0.990207 bsl -51. 

7916 

csl 0.993074 

esl 

0.0111668 



anl 

ooooo 

0 

10 0 0 

bnl 

ooooo 

2 

0 0 0 0 

as2 

0.431310 bs2 -59. 

6554 

cs2 3.503350 

es2 

0.0158715 



an 2 

0 10 0 0 

0 

0 0 0 0 

bn 2 

ooooo 

0 

0 2 0 0 

as 3 

0.794709 bs3 -113. 

2080 

CS3 3.168370 

es3 

0.0269699 



an3 

ooioo 

0 

0 0 0 0 

bn 3 

ooooo 

0 

0 0 2 0 

as4 

-0.652939 bs4 -9. 

,8232 

CS4 3.930330 

es4 

-0.0142865 



an 4 

oiooo 

0 

10 0 0 

bn 4 

ooooo 

0 

0 0 0 2 

as5 

1.158882 bs5 -76 

.8472 

cs5 8.532155 

es5 

0.0463471 



an5 

0 10 0 2 

0 

0 0 0 0 

bn 5 

ooooo 

0 

0 0 0 4 

as6 

0.980875 bs6 68 

.4453 

CS6 -10.5938 

es6 

-0.0414570 



an6 

oiooo 

0 

0 0 0 0 

bn6 

0 0 0 2 0 

0 

0 0 0 0 


0 2 
0.000 

0.000 


0 0 

0 2 

0 . 000 

0.000 


0 0 

0 2 

0.000 

0 . 000 


dsl -0. 

0 0 
0 0 
ds2 -0 

0 0 
0 0 
ds3 -0 

0 0 
0 0 
ds4 0 

0 0 
0 0 
ds5 -C 

0 0 
0 0 
ds6 C 

2 0 
0 0 


0.000 

0 . 000 

0.000 

1.000 

0.000 

0.000 

.343428 
. 340016 
.443814 
.163490 
1.868320 
). 574260 
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